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ABSTRACT 

We investigate the form of the one-point probabihty distribution function (pdf) for 
the density field of the interstellar medium using numerical simulations that successively 
reduce the number of physical processes included. Two-dimensional simulations of self- 
gravitating supersonic MHD turbulence, of supersonic self-gravitating hydrodynamic 
turbulence, and of decaying Burgers turbulence, produce in all cases filamentary density 
structures and evidence for a power-law density pdf with logarithmic slope around —1.7. 
This suggests that the functional form of the pdf and the general filamentary morphology 
are the signature of the nonlinear advection operator. 

These results do not support previous claims that the pdf is lognormal. A series of ID 
simulations of forced supersonic polytropic turbulence is used to resolve the discrepancy. 
They suggest that the pdf is lognormal only for effective polytropic indices 7 = 1 (or 
nearly lognormal for 7 7^ 1 if the Mach number is sufficiently small), while power laws 
develop for densities larger than the mean if 7 < 1. We evaluate the polytropic index for 
conditions relevant to the cool interstellar medium using published cooling functions and 
different heating sources, finding that a lognormal pdf may occur at densities between 
10^ and at least 10^ cm"^. 

Several applications are examined. First, we question a recent derivation of the IMF 
from the density pdf by Padoan, Nordlund & Jones because a) the pdf does not contain 
spatial information, and b) their derivation produces the most massive stars in the voids 
of the density distribution. Second, we illustrate how a distribution of ambient densities 
can alter the predicted form of the size distribution of expanding shells. Finally, a brief 
comparison is made with the density pdfs found in cosmological simulations. 



Subject headings: ISM: clouds - instabilities - magnetohydrodynamics - turbulence - 
STARS: IMF - cosmology 



1 



1. Introduction 

In order for stars or groups of stars to form, den- 
sities nmst be sufficiently large, and so the statisti- 
cal properties of the density field in the interstellar 
medium must be intimately connected with the star 
formation behavior of galaxies. The density field is 
coupled in a nonlinear manner to the velocity distri- 
bution of the gas, so the statistics of the density field 
can potentially serve as a diagnostic between differ- 
ent physical processes that might play a role in con- 
trolling the dynamics of interstellar material. For 
example, magnetically-dominated turbulence might 
exhibit significantly different density statistics than 
non-magnetic turbulence. Thus the density distri- 
bution could shed light on basic unsolved questions 
concerning star-forming regions, such as the nature 
and maintainance of highly supersonic motions ob- 
served at all but the smallest scales. The form of the 
density statistics might also be important in control- 
ling the ability of galactic gas to spontaneously de- 
velop hierarchical structure, as suggested by Vazquez- 
Semadeni (1994). Padoan, Jones & Nordlund (1997) 
have shown how observations of extinction fluctu- 
ations by Lada et al. (1994) can be used to con- 
strain the statistics of the density field, while Padoan 
(1995) and Padoan, Nordlund & Jones (1997) have 
attempted to relate the density distribution function 
to the stellar initial mass function, assuming a specific 
model for star formation. Understanding the physics 
behind the one-point density distribution function is 
also important in a cosmological context, since the 
form of this function in the nonlinear regime may be 
sensitive to the initial fluctuation statistics (see sec. 
4.4. below). 

In the present paper we focus on the behavior of 

the one-point probability distribution of gas densi- 
ties in numerical simulations relevant to galactic gas. 
This function, denoted /(p), is deflned here such that 
f{p)dp measures the fractional volume occupied by 
gas in the density range (p, p + dp), so f{p) is a prob- 
ability density function (pdf ) with units 1 / density. If 
the density pdf is parameterized in some range of den- 
sity as a power law with index 9, f{p) ~ p^ , then the 
fractional volume of space at density p per unit log- 
arithmic density interval and the fractional mass of 
material per unit density interval would be a power 
law with index (^ -|- 1). 

Vazquez-Scmadcni (1994) presented evidence from 
two-dimensional numerical simulations of randomly- 



forc:ed turbulence that the density pdf may be log- 
normal in form; i.e. that the logarithm of the density 
has a probability density that is normal (Gaussian). 
However, his simulations were purely hydrodynamic 
(i.e., neglected self-gravity, the magnetic fleld and 
the Coriolis force), isothermal, and were restricted 
to relatively small Mach numbers compared to what 
is observed in most interstellar regions. Vazquez- 
Semadeni, Passot fc Pouquet (1995) have presented 
simulations incorporating all of these physical ef- 
fects, but discussed the density pdf only in pass- 
ing. Nordlund & Padoan (quoted in Padoan et al. 
1997a, b,c) have also found a lognormal density pdf in 
randomly forced isothermal three-dimensional simu- 
lations. These simulations include the magnetic field, 
and reach higher Mach numbers, but contain no other 
astrophysically relevant physical ingredients. 

Part of the purpose of the present work is to ex- 
amine whether this lognormal density pdf persists in 
simulations that include most of the major physical 
processes expected to play some role in galactic turbu- 
lence, including self-gravity, stellar heating, the mag- 
netic field, rotation, and explicit heating and cooling 
functions in the energy equation, and which are car- 
ried out at large Mach numbers. If the results are 
different, we want to find out why, and try to iso- 
late the dominant physical effect (s). In fact, Gotoh & 
Kraichnan (1993) have presented numerical evidence 
(and an interpretation in terms of the "mapping" clo- 
sure) that the pdf for the Burgers equation (which is 
pressureless) is a power-law at high densities, rather 
than a lognormal. 

In section 2.1 we first show that simulations of the 
most complex system do not result in a lognormal 
pdf, but instead show evidence for power law behav- 
ior at densities above the mean. In order to under- 
stand the reasons for the different result and to isolate 
the physical process(es) that are most responsible for 
this behavior, we next successively reduce the system 
by removing physical processes. In sec. 2.2 we exam- 
ine simulations similar to the original set, but without 
magnetic fields and rotation. Again we find the power 
law pdf, with a similar slope and density range as be- 
fore. As a more radical reduction of the system, in sec. 
2.3 we study the density pdf for simulations that con- 
tain no physics at all except for nonlinear advection in 
the continuity and momentum equations (Chappell & 
Scalo 1997). These simulations contain no pressure, 
and so are equivalent to a gas with effective polytropic 
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index equal to zeroj^ Surprisingly, the spatial density 
fields and the density pdfs resemble the earlier simu- 
lations at densities above the mean, suggesting that 
advection is the dominant process. The importance 
of the advection term is well known in incompressible 
turbulence, in which even the pressure term can be 
incorporated into the nonlinear operator. 

In order to understand these results, in section 3 
we present a few conclusions that can be drawn from 
a large study of one-dimensional forced turbulence 
by Passot & Vazquez-Semadeni (1997), who have in- 
vestigated the statistics of the density and velocity 
fields for a wide range of Mach numbers and poly- 
tropic indices. It turns out that a lognormal density 
pdf should only obtain at small Mach numbers or, for 
any value of the Mach number, when the polytropic 
index is equal to unity. This explains the discrepancy 
between the present simulations (all of which have 
small effective polytropic indices) and previous work 
(which assumed 7 = 1). We also discuss the expected 
value of the polytropic index as a function of density 
and temperature, based on published cooling func- 
tions and various heating sources and conclude that 
at densities less than about 10'^ cm~'^ the polytropic 
index should be significantly less than unity, favor- 
ing power-law density pdfs. For higher densities we 
show that saturation of the cooling function should 
lead to 7 « 1 and therefore possibly a lognormal pdf. 
However at still larger densities 7 gas-grain collisional 
cooling or heating should dominate, forcing 7 away 
from unity again. Some implications of these results 
are discussed in section ^ In particular, we argue that 
the distribution of masses of collapsing structures or 
stars (the IMF) cannot be derived from the one-point 
density pdf, and illustrate how a distribution of den- 
sities might affect the predicted size distribution of 
expanding shells driven by young stars. 

2. Three types of simulations and their pdfs 

2.1. Self-gravitating magnetohydrodynamic 
models with rotation 

A first class of models considered in the present 
paper consists of the fully nonlinear self-gravitating 
magnetohydrodynamic (MHD) equations with added 



^The polytropic index characterizes the response of the gas to 
quasi-static compression or expansion in thermal equilibrium, 
and is of course different from the ratio of specific heats which 
is, for example, 5/3 for a monotonic gas with no internal degrees 
of freedom. 



model source terms for radiative cooling, stellar and 
diffuse heating, and the Coriolis force, as described 
in Passot, Vazquez-Semadeni & Pouquet (1995). The 
equations of this model, to which we will refer as the 
MHD model, are: 
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A = p'L{T), (9) 



LiT) = L,T'^ for T,<T< T,+i, (10) 



Li = 1.14x10^5 bi = 2 

L2 = 5.08 X 10^6 62 = 1.5 

L3 = 2.35 X 10" 63 = 2.867 

La = 9.03 X 10^8 64 = -0.65 



where 

Ti = 100 
r2 = 2000 
T3 = 8000 
Ti = 10^ 
T5 = 4 X 10^ 

We refer the reader to Passot et al. (1995) and to 
Vazquez-Semadeni, Passot & Pouquet (1996) for a 
detailed discussion of this model. Here we just de- 
scribe it briefly. The numerical method used for solv- 
ing these equations, as well as the purely hydrody- 
namic (HD) runs discussed in i |2.2| , is pseudospectral, 
implying that diffusion operators have to be included 
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explicitly, since the method does not produce any nu- 
merical viscosity. The initial conditions are Gaus- 
sian with random phases and a spectrum of the form 
P{k) oc fc"* exp(— fc^/16) in all variables. All the simu- 
lations presented in this section are two-dimensional 
(2D). 

In the momentum equation, eq. (j^), f2 is the an- 
gular velocity of Galactic rotation, appearing in the 
Coriolis term, and vs is a hyperviscosity coefficient 
also used in the magnetic flux freezing equation, eq. 
(|^). The usage of hyperviscosity of the form V* in- 
stead of a standard Laplacian viscosity operator con- 
fines viscous effects to the very smallest scales in the 
simulations, allowing the development of larger tur- 
bulent inertial ranges. However, it produces oscilla- 
tions in the vicinity of strong shocks, as discussed in 
S3.1, and thus a small amount of second-order vis- 



cosity has been added in some runs. A mass diffusion 
term, with a corresponding coefficient is used in the 
continuity equation, in order to smooth out the den- 
sity gradients, thus allowing the simulations to reach 
higher r.m.s. Mach numbers. A table with the full list 
of fiducial parameter values is presented in Passot et 
al. (1995). 

In the internal energy equation, the terms Fd and 
Fg respectively refer to diffuse background and stel- 
lar heating rates per unit volume. The latter mimics 
the heating by ionization occuring in the HII regions 
surrounding OB stars, which are modeled as point 
heating sources appearing at every location where the 
density exceeds a threshold which we arbitrarily set 
at a density pc ~ 30 (p). The "stars" have a lifetime 
of 6.5 X 10^ yr. Finally, A parameterizes the radiative 
cooling rate per unit volume. Note that, for com- 
patibility with standard notation, we have changed 
the notation from that in Passot et al. (1995) and 
Vazquez-Semadeni et al. (1996). In those papers, Fs, 
Fd and pA were the heating and cooling rates per unit 
mass, and the exponent in the equation defining Fd 
was labeled a, being related to the exponent used in 
the present paper by a = 1 ~ a. 

Figure ^ shows a typical view of the density field 
for one MHD two-dimensional (2D) simulation at a 
resolution of 800 x 800 grid points at time t = 1.2io, 
where to = 8.2 x 10^ yr is the sound crossing time 
for the integration domain. In this run, the dissipa- 
tion coefficients are — 5.63 x 10^^^, = = 0, 
p = 3.28 X 10"'^, and a = 1/2. The box size corre- 
sponds to a region of size 1 kpc in the Galactic plane 
at the solar galactocentric distance. We refer to this 



simulation as MHD800. It is actually a restart of the 
run named "run 28" in Passot et al. (1995), but at 
larger resolution (run 28 had 512 x 512 grid points). 
Moreover, the star formation has been turned off in 
the last stages of this simulation in order to allow the 
turbulence to develop larger density peaks. Other- 
wise, the star formation criterion, which turns on a 
star whenever the density exceeds pc, tends to prevent 
the simulation from reaching densities significantly 
higher than pc, since the stars heat their surround- 
ings, increasing the pressure and producing expand- 
ing bubbles, and reducing the density. This run is 
discussed in Vazquez-Semadeni, Ballesteros-Paredes 
& Rodriguez (1997), where another view of the den- 
sity field att = 0.9to is presented. Note that the only 
stellar energy injection considered is that due to ion- 
ization heating of the medium surrounding OB stars, 
supernova explosions not being included. The latter 
implies that the largest complexes (see below) cannot 
be blown apart by the "stars" . 

The density field is seen to consist of dense, fila- 
mentary, knotty structures ( "clouds" ) which are inter- 
connected in extremely complicated patterns ("com- 
plexes"), being the local peaks within less dense, ex- 
tended yet highly amorphous larger structures ("dif- 
fuse clouds") (e.g., the two large structures in the 
upper- and lower- right parts of the figure). In turn, 
the diffuse clouds gradually disappear into the lowest- 
density "intercloud medium" (dark regions in the 
middle and left upper and lower parts of the image). 

Figures |^a and respectively show the histograms 
of logp {— pf{p), where / is the density pdf), for run 
MHD800 at two different times, h = OMo and t2 = 
1.2^0, the latter corresponding to the image shown 
in fig. |l|. The density is shown in units of the mean 
density in the field. The histogram is computed over 
the whole field, so it contains 800^ sample points. 

The histogram at ti exhibits a clear power-law be- 
havior in the interval < logp < 1.3, with slope 
—0.73. The histogram at t2, on the other hand, 
does not exhibit such a clear power-law behavior, al- 
though the region < logp < 1.5 can be roughly 
approximated by a power-law of slope —0.63. Both 
power-laws are indicated by the solid lines. Note that 
these slopes correspond to density pdfs / with slopes 
steeper by one power of p, i.e., / oc p^^ '^. 

For densities outside the ranges mentioned above, 
both histograms turn over, decaying rapidly at both 
very large and very small densities. The drop-off at 
large densities can be easily understood as a conse- 
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quence of the viscosity and mass diffusion included in 
the equations (the and /i terms in eqs. and j^], 
respectively), whose effect is to smooth out tall, nar- 
row density peaks, "losing" them from the histogram. 
We speculate that the self-similar range would extend 
further into higher density values if more resolution 
(and consequently lower dissipative coefficients) could 
be used. 

The turnover at low densities occurs close to the 
mean density, but has a much more abrupt charac- 
ter than either the Burgers '-type runs discussed in § 
2.2 or the one-dimensional (ID) simulations of Passot 
& Vazquez-Semadeni, discussed briefly here in § p?l] . 
This is probably due to the presence of self-gravity in 
run MHD800, which has managed to collect the gas 
in large complexes, from which the gas cannot escape. 
Therefore, the gas in the voids is only slowly gravi- 
tationally accreted into the complexes, but cannot be 
swept up by the filaments as in the Burgers and ID 
runs. That is, the density minima in run MHD800 are 
probably of gravitational, not turbulent origin, thus 
taking much longer times to be evacuated. 

2.2. Hydrodynamic self-gravitating runs 

In this section we briefly discuss a first simplifica- 
tion step down from the full eqs. (||) to (||), obtained 
by eliminating the magnetic fleld and the Coriolis 
force due to Galactic rotation in the simulations. This 
is an important case, because these physical agents 
can be the source of added "hardness" in the flow 



(§|3.2|; see also Vazquez-Semadeni et al. 1996). 

Figures ^a and |^b show the density field for a 2D 
run (labeled HD512) similar to run MHD800 except 
without magnetic fields or the Coriolis force, and at 
a resolution of 512 grid points per dimension. At this 
lower resolution, the values of the diffusion parame- 
ters are i/g = 10"^^, 1^2 = 2 x lO'^, /i = 8 x 10"^, and 
a — 1/2. The run considered evolves in time as fol- 
lows. The initial turbulent transients induce the for- 
mation of clumps and filaments, which over the first 
0.48 crossing times of the simulations (= 3.9 x 10^ yr) 
merge into larger and denser structures due to the 
combined action of turbulence and self-gravity but 
with no star formation activity. After this period, 
star formation begins, producing isolated expanding 
"HII regions" . However, since no supernovae are in- 
cluded, the stellar energy input is insufficient to halt 
collapse of the large structures, and by t = 1.55to 
the largest structure finally collapses. Note that this 
type of collapse did not occur in the simulations of 



Vazquez-Semadeni, Passot & Pouquet (1995), which 
were also non-rotating and non-magnetic, because the 
density threshold for star formation used in that pa- 
per was much lower (8(p)) and the fluid was harder 
(a=l). 

The density fleld is shown at i = OAto (fig. ||a) and 
t = 1.47<o (fig- ^)- The first time corresponds to 
an epoch shortly after the initial transients, at which 
matter has not significantly gathered gravitationally, 
and star formation events have not yet begun. The 
corresponding histogram for log p for this epoch is 
shown in fig. and contains 512^^ data points. Again 
a near power-law with slope ^ —0.6 can be seen in 
the range — 0.7 < logp < 0.4. Interestingly, at larger 
densities (0.4 < logp < 1.5) another power-law can 
be fitted, with slope ^ —2.0. Thus, power-law ranges 
in the density pdf appear to be present even after the 
removal of magnetic fields and rotation. 

The density field shown in fig. |^b corresponds to 
the final state of run HD512, in which a large struc- 
ture has collapsed gravitationally in the upper left 
quadrant of the field. The maximum density for 
this field, occuring at the center of the collapse, is 
p — 703 (p). The corresponding histogram for logp 
is shown in fig. ^jb. Quite surprisingly, even for this 
highly singular density configuration, the logp his- 
togram exhibits a power-law range, although with 
two significant bumps at 1 < logp < 1.5 and at 
2.5 < logp < 2.75. These bumps thus seem to corre- 
spond to two special densities: the density threshold 
for star formation (recall pc = 30 (p)), and the density 
of the collapsed region. The bump at pc can be eas- 
ily understood since the stellar heating tends to pre- 
vent the density from exceeding pc, except when the 
self-gravity of the complex is so large that the stellar 
heating is insufficient to prevent generalized collapse. 
Note again that this would not occur if supernovae 
capable of blowing away the complex were included. 
Nevertheless, the rest of the histogram continues to 
be a power-law. The measured slope is —1. Thus, the 
power-law behavior at large densities is preserved in 
the absence of magnetic fields and rotation, and even 
in the presence of local gravitational collapse events. 

2.3. Stripped-down quasi-Burgers simulations 

A third series of simulations was run with the in- 
tention of stripping-down the momentum equation to 
isolate the effects of the advection term u ■ Vu in 
the momentum equation. For this purpose the equa- 
tions of continuity and momentum conservation were 
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solved, but all terms on the RHS of the momentum 
equation were set to zero, except for an additional 
stellar forcing term equivalent to that used in the 
MHD and HD simulations. Both forced and decay- 
ing simulations were performed. Thus these calcu- 
lations describe a highly compressible fluid in which 
advection and the corresponding "ram pressure" com- 
pletely dominate the thermal pressure, or, equiva- 
lently, in which the effective polytropic index 7 is zero. 
The system is equivalent to a forced Burgers flow, ex- 
cept that the viscosity in our case is numerical diffu- 
sion and the forcing is coupled spatially to the density 
field, since star formation and momentum input are 
assumed to occur at a theshold density. However, the 
runs with threshold star formation were not capable 
of generating a sufficiently broad range of densities 
at the high density end to learn much about the pdf, 
so we forego any discussion of the forced cases and 
consider only the pure decay simulations. The runs 
with star formation (most of the runs) will be dis- 
cussed elsewhere (Chappell & Scalo 1997). Because 
of the similarity to Burgers turbulence, we refer to 
these simulations as series "B," or "quasi-Burgers" 
runs. 

The advection terms were differenced according to 
a variant of a Van Leer (1977) first-order scheme, al- 
tered to eliminate most of the spurious anisotropy of 
numerical diffusion in that scheme. The boundary 
conditions were doubly-periodic. The initial condi- 
tions were a uniform density field and a Gaussian ve- 
locity field with prescribed power spectrum. Runs 
with resolution of 128^, 256^, and 512^ were exam- 
ined. The scales were normalized such that the lattice 
spacing was 7.8 pc, so these resolutions correspond to 
a total region size of 1, 2, and 4 kpc, although the 
adopted size normalization is probably irrelevant for 
the results presented here concerning the density pdf. 
For pure advection, the hydrodynamic equations are 
invariant with respect to a change of scale, so these 
simulations can be applied to any scale. 

All the simulations develop into a network of irreg- 
ularly shaped filaments which cover a range of sizes. 
These filaments are the products of the advection 
operator nonlinearly self-organizing the initially ran- 
domized velocity field into a collection of shocks. The 
thicknesses of the filaments are set by numerical dif- 
fusion; they are thinner than in the MHD and HD 
cases because there is no pressure force. It is just this 
absence of pressure that ensures the larger compress- 
ibilities that feed the filament-generating advection 



operator. With no momentum input, the filaments 
simply continue to sweep up material and each other, 
so that as time proceeds the structure becomes more 
concentrated on large scales and the velocities mono- 
tonically decrease. An image of the density fleld (ac- 
tually of p^/^) at two times is shown in Fig. ^. It is 
seen that the structure resembles the HD simulation 
shown in Fig. |^a at a time early enough that self- 
gravity has not yet caused the fllaments to become 
concentrated into "complexes." 

Figure ^ shows the development of the density pdf 
(actually pf{p)) for two sets of 256'^ simulations, for 
initial conditions with velocity energy spectra propor- 
tional to fc~^ (left) and fc" (right). Each pdf is based 
on about 10^ points, but most of them are in the 
voids between the filaments. The solid line is a ref- 
erence lognormal pdf with peak at p — 1, shown just 
to illustrate the degree to which the simulation pdfs 
do or do not match a lognormal. Times are, from top 
to bottom, in units of initial crossing times, 0.006, 
0.07, 1.1, and 17. After about 1 crossing time one 
can see the same sort of power law for log p ^ that 
was found in the MHD and HD simulations. For both 
sets of initial conditions the logarithmic slope in the 
density is about —0.7. The pdfs become noisy at late 
times because the density field has evolved into a rela- 
tively small number of large scale filaments. Although 
the number of points is too small to reveal the nature 
of the highest-density tail, the dropoff seen in the k~'^ 
case, and to a lesser extent in the fc*^ case, is consistent 
with the behavior seen in the better-sampled MHD 
and HD pdfs; the dropoff at highest densities is once 
again probably due to viscosity, which smears out 
the smallest scale structures, which are the thinnest, 
densest filaments that have recently undergone a col- 
lision with another filament. 

We place no physical significance on the low den- 
sity (p ^ 1) part of the pdf. In this pseudo-Burgers 
flow fllaments sweep up all the inter-fllament material, 
and once the entire simulation area has been crossed 
by at least one fllament the density between fllaments 
would be zero except for the action of viscosity (nu- 
merical diffusion here) which causes material to leak 
back into the "voids." The low-density pdf extends to 
much lower densities than shown in the figure. 

Figure shows the late-time pdfs for three 256^ 
simulations initially excited at a single given wavenum- 
ber fcp = 64, 16, or 4. Although the impression is 
subjective, it appears that the pdfs are better rep- 
resented, for (0 2, by a power law plus a steeper 
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dropoff at large densities. The "knee" that sepa- 
rates the power-law density range from the dropoff 
appears to shift towards lower densities with increas- 
ing ko, occurring at p « 60 for ko — A and p « 30 
for ko = 16. In both cases the logarithmic slope of 
the power law is around -0.6 to -0.8. For ko = 64 
the pdf appears to be fit by the lognormal for p^2, 
but inspection of other plots for this case at differ- 
ent times suggests that the knee has moved to such 
small densities, p ~ 10, that the power law region 
above p ~ 2 spans too small a range to be clearly 
visible. The knee in the ko — 64 case becomes more 
prominent at times later than shown in Fig. This 
is because, when the initial velocity fluctuations have 
such a small scale, it takes many filament interactions 
to build up large structures. The same effect can be 
seen in Fig. ^ where the power law develops much ear- 
lier for the k~'^ initial energy spectrum than for the 
spectrum. 

In order to increase the sampling of cells at large 
densities, we ran a 512^ simulation with a broad-band 
initial spectrum. The resulting density pdf is shown in 
Fig. ^ at two times. At early times the pdf at densities 
greater than the average appears to fit a lognormal, 
but after about 1.5 crossing times the result is more 
suggestive of a short power law segment with a falloff 
at p ^ 10. Between logp = 0.2 and 1.0 the logarithmic 
slope is about -0.6. 

We note that Gotoh and Kraichnan's (1993) 1- 
dimensional simulations of the decaying Burgers equa- 
tion produced a density pdf with a clear power-law 
form, although steeper than found here. Their power 
law regime had a larger extent in density probably be- 
cause of the much greater dynamical range possible in 
ID compared to the present 2D simulations. We take 
their work as support for the power law pdf regimes 
suggested here. 

Because of the poor sampling statistics for these 
quasi-Burgers decay simulations compared to the MHD 
and HD runs, the results are not conclusive, but the 
pdfs are generally consistent with the power law seg- 
ment plus high density falloff found in the more phys- 
ically complete MHD and HD simulations. The rough 
agreement is especially surprising because the B sim- 
ulations used a very different numerical method and, 
in most cases, different initial conditions. Since the 
B simulations have omitted all the physics except for 
nonlinear advection, this rather surprising agreement 
suggests that the physical process responsible for the 
power law portion of the pdf in the high- density side 



is nonlinear advection. Intuitively, this can be under- 
stood in terms of the tendency of the advection term 
to induce multiplicative processes which may lead to 
power laws, if the equation of state permits it (Passot 

6 Vazquez-Semadeni 1997, hereafter PVS). 

3. The Role of the Polytropic Index 
3.1. One-dimensional simulations 

The simulations discussed in the previous sections 
suggest (with varying degrees of precision) the devel- 
opment of power-law regions at high densities in the 
density pdf of interstellar gas. This result is in con- 
trast with previous work (Vazquez-Semadeni 1994; 
Padoan et al. 1997a, b,c) finding lognormal distribu- 
tions in a variety of numerical simulations. Vazquez- 
Semadeni (1994) presented isothermal (7 = 1) sim- 
ulations of two-dimensional, randomly-forced Navier- 
Stokes turbulence in the weakly compressible regime. 
No other processes, such as star formation, self- 
gravity, magnetic fields, etc., were included. Padoan 
et al. (1997a, b, c) refer to highly-compressible three- 
dimensional isothermal MHD simulations with ran- 
dom forcing and without self-gravity to be presented 
elsewhere, although, to our knowledge, the pdfs have 
not been published yet. For both sets of simulations, 
a lognormal pdf is reported. 

The discrepancy between our pdfs and those of 
Vazquez-Semadeni (1994) and Padoan et al. (1997a, b,c) 
can be understood in terms of some recent results by 
PVS, which we briefly summarize here. As a first 
step towards understanding the relation between the 
dynamical and statistical properties of highly com- 
pressible turbulent flows, PVS have performed a sys- 
tematic investigation of one-dimensional (ID) simu- 
lations of pure Navier-Stokes polytropic turbulence 
(with mass diffusion added on some runs in order to 
allow the simulations to tolerate strong shocks), find- 
ing that a lognormal density pdf only appears when 

7 = 1. In fact, as explained below, PVS suggest that 
the case 7 = 1 may be singular in this respect. 

The equations solved by PVS in nondimcnsional 
form are: 

dp d{pu) d'^p 
dt dx 

and 



(11) 



du du 
dt dx 



1 



+ (12) 



jM'^p dx dx 
where M is the Mach number of the velocity unit 
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In these simulations, only second-order viscosity is 
used in order to avoid the spurious oscillations in- 
duced in the vicinity of shocks by the use of hypervis- 
cosity (Passot & Pouquet 1988). Mass diffusion (the 
RHS of eq. ) is used in small amounts to help the 
simulations survive strong shocks. The probability 
distributions of log p are obtained by considering all 
grid points (typically 2048) and averaging over nearly 
150 code non-dimensional time units, sampling every 
0.1 time units. Thus, the histograms contain over 3 
million points. 

As an illustration, we show in figs. and b the 
density fields of two different simulations, both with 
M = 3, and 1^2 = A* = 0.003, but with 7 = 1 in (a) 
and 7 = 0.3 in (b). These runs start at rest, but are 
driven with a random acceleration at wavenumbers 
2-20, with a correlation time <corr = 4.8 x \0~Ho- The 
associated histograms for log p are shown in figs. |l^ 
and b. The dotted lines show a least-squares fit to 
a lognormal curve. As can be seen, for the simula- 
tion with 7 = 1 the fit is excellent, suggesting that 
indeed the logp histogram is a true lognormal. In- 
stead, in figure |lo|b, a clear power-law tail is seen to 
develop at large density fluctuations. Furthermore, in 
fig. |ll] we show the corresponding histogram for a run 
with 7 — 0.3 but AI = 1.8. Upon comparison with 
fig. ^p, this run illustrates the effect of varying the 
Mach number M at a fixed 7. The deviation from 
a lognormal towards a high-density power-law tail is 
more noticeable at large M . These results are con- 
sistent with the power-law functional form of the pdf 
at high densities found by Gotoh & Kraichnan (1993) 
in their ID numerical simulations of decaying Burg- 
ers turbulence, which effectively have 7 = 0. Also, 
note that the development of the high-density power- 
law tail in the pdfs of the present simulations cannot 
be attributed to a numerical effect, such as the in- 
clusion of the mass diffusion term in the continuity 
equation. Runs at lower Mach numbers (not shown), 
which could be performed without the need for such a 
term, still exhibited a power-law tail, although not as 
well developed, due to the Mach number dependence 
discussed above. 

In order to interpret these results, it is instruc- 
tive to rewrite the inviscid one-dimensional gas dy- 
namics equations in Riemann invariant form. The 
Riemmann invariants of the problem arc defined, for 
7 7^ 1, as = u:kKc^ where u is the one-dimensional 
fluid velocity, c = p'^^^^^l'^ jM is the sound speed, and 
K = ±2/(7 — 1). They are just advected along the 



characteristics whose speeds are u ± c. In the case 
J = 1, — u ± log p/M while the characteristic 
speeds are u ± l/M. One can easily see that, under 
the transformation p^l/p, 7— i'2 — 7, the charac- 
teristic speeds remain unchanged while the Riemann 
invariants are exchanged. Under such a transforma- 
tion the local dynamics is not preserved but the sta- 
tistical quantities such as the pdf remain close. For 
7=1 this indication is consistent with the fact that 
the pdf is a lognormal, a curve symmetric under the 
change log p — log p. 

It should be noted that the lognormal pdf at 7 = 1 
disappears when a random force fr is used, leading 
to replace by fr/p in eq. (p^). We speculate that 
this is due to the strong nonlinearity introduced by 
the division by p. Thus, to study the effect of the 
advection operator on the density pdf, the usage of 
an acceleration seems more adequate. 

In terms of these results, we can interpret the pre- 
vious claims of lognormal pdfs (Vazquez-Semadeni 
1994; Padoan et al. 1997a, b,c) in a simple way. Since 
both sets of simulations were isothermal, their result- 
ing pdfs were bound to be lognormal, independently 
of the Mach number. However, for any other values 
of the polytropic exponent, power-law pdfs should be 
realized. 

3.2. The polytropic index of cool interstellar 

gas 

The one-dimensional calculations of supersonic 
Navier-Stokes turbulence discussed above and pre- 
sented in detail by PVS strongly suggest that the form 
of the density pdf f{p) depends sensitively on the ef- 
fective polytropic index 7. In particular, if the turbu- 
lence is highly supersonic, f{p) is strictly lognormal 
only for 7 equal to unity, while for smaller (larger) 7 a 
power law develops for densities larger (smaller) than 
the mean. The sensitivity of hydrodynamic behavior 
to 7 has been found previously in non-turbulent con- 
texts. For example, Pongracic (1994) and Foster & 
Boss (1996) showed that the ability of a shock inci- 
dent on a model cloud to induce gravitational collapse 
depends on 7 (through a piecewise power law cooling 
rate for Pangracic). Vazquez-Semadeni et al. (1996) 
estimated a critical value of the polytropic exponent 
below which turbulence can induce bound condensa- 
tions, as a function of the dimensionality of the tur- 
bulent compressions. 

The polytropic index 7 defined by P ~ p''', is iden- 
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tically zero for the quasi-Burgers simulations. For hy- 
drodynamic simulations which allow for heating and 
cooling, 7 can be determined locally by the tempera- 
ture and density dependence of the heating and cool- 
ing function. If the heating and cooling rates per unit 
volume depend on temperature and density according 
to r p"^ and A ^ p'^T^, then the condition of ther- 
mal equilibrium (which holds for the flows of interest 
because of the small cooling times, except for just 
behind shocks), gives 

^ ^ ^ ^ dlogT ^ ^ _ jc-a) / _ logp db \ 
d log p b \ b d log p J 

(13) 

(cf. Elmegreen 1991; Passot et al. 1995; Vazquez- 
Semadeni et al. 1996). The last factor in parentheses 
allows for the possibility that the temperature depen- 
dence exponent b in the cooling function is itself a 
continuous function of the density (see eq. ||l^ be- 
low). 

In the MHD and HD simulations discussed above, 
the heating was assumed to be due to a combination 
of stellar and diffuse radiative heating, with a = 0.5 
for the latter, while the cooling rate was a piece- 
wise power-law function appropriate to temperatures 
larger than a few hundred degrees, because the sim- 
ulations were designed to model large-scale dynamics 
of the atomic gas. For those simulations the values 
of 7 are in the range 0-0. 5. Q Therefore the develop- 
ment of the power law f{p) at large densities found in 
the MHD, HD, and B simulations is consistent with 
the one-dimensional results, since 7 is significantly 
smaller than unity. 

However, the most interesting applications of the 
density pdf concern regions that can form stars, which 
are at small temperatures and are generally molecu- 
lar instead of atomic. We therefore need to estimate 
the value of 7 in cool molecular gas. A few previous 
attempts in this direction are summarized by Larson 
(1985, Fig. 2) who suggests 7 « 0.7 in the absence of 
stellar heating, for T ~ 10 - 50 K and n ~ 10^ - 10^^ 
cm~'^. However there are a number of possible heat- 
ing mechanisms, and the cooling is complicated by the 
contribution of different coolants and optical depth 
effects at large densities. 

A major study of interstellar cooling and heating 
at low temperatures (< lOOK) was given by Gold- 

Actually, the range of cooling functions described in sec. ^ give 
7 = 3.3 for 10® < T < 4 X 10^ K. However, those temperatures 
are never reached in the simulations. 



smith and Langer (1978, hereafter GL) for the density 
range 10^-10^cm~'^. The resulting cooling function is 
dependent on the prescribed abundances of coolants 
and the treatment of radiative transfer (escape prob- 
ability method with constant large velocity gradient) . 
A considerably more detailed treatment of the cooling 
function was presented by Neufeld, Lepp, and Melnick 
(1995, hereafter NLM) for densities lO^-lO^Ocm-^ 
and temperatures 10-2500K. Besides including more 
coolants, using updated data on cross sections, and 
a modified velocity gradient parameter for the radia- 
tive transfer, NLM obtained coolant abundances by 
a self-consistent (but steady-state) solution to a large 
chemical reaction network. The results assume com- 
plete shielding from ultraviolet radiation, and may be 
significantly different in lower-column density regions. 
We base the following discussion of 7 on the results 
of GL and NLM, concentrating on densities 10^-10^ 
cm~'^ and temperatures 10-100 K. We emphasize that 
our purpose is not to give a comprehensive treatment 
of the cooling and heating functions, but only to out- 
line what seem to be the dominant physical effects on 
7 and to motivate a more detailed treatment. 

For temperatures less than about 100 K, we find 
that the temperature dependence of the cooling func- 
tion, parameterized by L ~ T^, can be adequately 
represented by a density-dependent function, based 
on Table 4 of GL, 

6=i(l + logn) (14) 

for log n in the range 2 to 5, where n is the total parti- 
cle density. The numerical and graphical results pre- 
sented by NLM are not sufficient to confirm this rela- 
tion, although the trend of larger temperature depen- 
dence at larger densities seems to hold for densities 
up to 10^ cm^^, based on Fig.2 in NLM. At smaller 
densities, based on the data in GL, b is smaller than 
given by eq. (0). For example, at n — 10^ cm~^, 
b « 0.5. 

For densities less than a few hundred cm~'^ the to- 
tal cooling rate per unit volume varies with the square 
of the density (c = 2), but at densities around 10'^ 
cm^'^ or larger, the density dependence becomes much 
weaker because of radiative trapping (e.g. GL). The 
density at which the flattened density dependence oc- 
curs depends on a number of effects, but is about 
lO'^cm"'^ for the parameters adopted by GL. The re- 
sulting density dependence is uncertain because of the 
approximate nature of the radiative transfer calcula- 
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tions in the available work and the assumption of a 
smooth (not turbulent) velocity field. Inspection of 
Fig. 3 of NLM indicates that for n between 10'^ and 
10^ cm -3, c is about 1.6 for T = 100 K and 1.2 for 
T = 10 K, although at T = 10 K c is very close to 1.0 
for their singular isothermal sphere model. For rea- 
sons we do not completely understand, the density 
dependence exponents at high densities are smaller 
by about 0.6 in the calculations of GL, Figs. 7a-c, 
with c as small as 0.5 at T = 10 K. Some possible 
reasons for the differences between GL and NLM are 
given by NLM, although they do not explicitly discuss 
the differences in density dependence. For illustrative 
purposes we will adopt c = 1.2 for n > 10'^ cm""^ and 
T = 10 K. 

The heating rate is problematic because it is un- 
certain what the dominant process is (see GL, Black 
1987, and HoUenbach 1988 for reviews). It is often 
assumed that the main heating agent is cosmic rays, 
with r ~ p, so that a = 1. The heating due to II2 
formation on grains depends on the fraction of neu- 
tral hydrogen, but at large densities the heating rate 
is again proportional to p. Heating rates due to gas- 
grain collisions, compression or collapse, ambipolar 
diffusion, or turbulence are more complicated and un- 
certain, and may include a temperature dependence. 
We ignore these sources here, except to note that b 
may be significantly larger than unity for these cases. 
An extreme case is gas-grain heating, which is pro- 
portional to p'^T^I'^{Tgr - T), where Tgr is the gram 
temperature. This case is discussed below. 

Our results suggest different behavior of 7 in four 
density ranges. 

1. n < 10"^ cm^"^. In this case, taking the density 
dependence of the heating as a = 1, and including the 
density dependence of 6 (eq. Ill]) in eq. (13), an effect 



that pushes 7 closer to unity, the polytropic index is 
0.56 for n = 10^ (6 = 1.5) and 0.88 for n = 10^ (6 = 
2). These bracket the estimate by Larson (1985). At 
densities as low as 10 cm~^, b « 0.5 based on the 
GL cooling rates, and 7 is negative, corresponding to 
thermal instability. 

We conclude that, as long as the dominant heat- 
ing mechanism is cosmic rays, H2 formation, or some 
other process whose rate scales with the density in a 
manner not much steeper than linearly, the polytropic 
index will be sufficiently smaller than unity in molec- 
ular regions with n < 10'^ cm~^ that the power law 
density pdfs found in the simulations should apply 
there, also. A lognormal pdf should only obtain for 



a nearly quadratic density dependence of the heating 
rate per unit volume. 

2. n ~ 10^ to a few times 10"' cm^"^. The situa- 
tion is different at densities larger than 10'^ cm~'^, at 
least when gas-grain collisional cooling does not dom- 
inate, because radiative trapping weakens the den- 
sity dependence of the cooling rate, i.e. reduces the 
exponent c. At T = 10 K, adopting c = 1.2 from 
the results of Neufeld et al. (1995) and again using 
a = 1 for the density dependence of the heating rate, 
7 varies from 0.98 at n = 10^ to 0.99 at n = 10^. At 
T = 100 K, taking c = 1.6, 7 varies from 0.92 to 0.97 
over this density range. The polytropic indices are 
close enough to unity that the lognormal density pdf 
might occur. There is some inconsistency in arriving 
at these numbers, because the density dependence of 
the temperature exponent in the cooling rate (eq. 14) 
was based on GL, while the density exponent is based 
on NLM. Without the density dependence of 7, the 
derived values of 7 are slightly smaller, but still close 
to unity. 

That 7 may be very close to unity at densities 
above 10^ cm~'^ is by no means a definitive result. Af- 
ter the above calculations were done, we found that 
Lis & Goldsmith (1990) have given polynomial fits 
to the cooling function of GL, which includes radia- 
tive trapping. Taking analytical derivatives of their 
expressions in order to evaluate 7, and assuming a 
linear dependence of the heating rate on density, we 
find that at log n = (2,3,4) 7 = (0.89, 1.21, 1.26). 
This suggests that 7 will only be very close to unity 
in a narrow density range around 300 cm""^. The rea- 
son for the difference is that the effective values of c, 
the density dependence of the cooling rate, are con- 
tinuously varying and significantly smaller than the 
value c = 1.2 which we estimated for the high-density 
cooling functions of NLM. 

Based on the NLM cooling rate, we suspect that a 
major transition occurs at densities above which ra- 
diative trapping in the cooling lines becomes impor- 
tant. At smaller densities the polytropic index should 
be significantly smaller than unity, and a power law 
density pdf should occur. At larger densities the value 
of 7 should be close to unity, and a lognormal density 
pdf is possible. Since these different density pdfs are a 
reflection of the velocity field (which generates density 
fluctations through the dilatation V • u), we expect 
that the velocity field and other associated phenom- 
ena should also be qualitatively different in the two 
density regimes. One caveat to this conclusion is that 
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radiative trapping may be very different for a turbu- 
lent velocity field than for the linear velocity gradients 
assumed in the existing cooling rate calculations. (For 
a calculation of molecular line formation in a turbu- 
lent velocity field, but without density fluctuations, 
see Kegel, Piehler & Albrecht 1993.) Another caveat 
is the fact that we obtain values of 7 significantly 
larger than unity for n = 10^-10''cm~^ using the fits 
of Lis & Goldsmith (1990) to the cooling function of 
GL, as discussed above. 

The possibility that 7 may be close to unity at den- 
sities between 10^ and 10** cm~'^ is consistent with 
the near constant temperature ~ 10 K observed in 
dark clouds without internal protostellar heat sources. 
However, this latter result is based largely on CO ob- 
servations that sample a fairly narrow range of den- 
sities. Often this result is ascribed either to the efh- 
ciency or the temperature dependence of the cooling; 
however it is clear that near isothermality should oc- 
cur in this density range primarily because the density 
dependence of the cooling rate changes as radiative 
trapping becomes important. 

3. n ^ 10^ cm~'^. At still larger densities, a differ- 
ent physical effect becomes important. If there are no 
embedded protostellar sources to heat the dust grains, 
then cooling by gas-grain collisions should dominate 
the molecular cooling at densities above 1 x 10** cm~^ 
to 5 X 10"* cm~'^, depending on the grain tempera- 
ture and grain parameters. (For the gas-grain energy 
transfer rate, see, for example, GL and HoUenbach 
and McKee 1989.) The gas-grain cooling rate is pro- 
portional to p'^T^/^{T - Tgr), so c = 2 and 6 w 3/2 
if T » Tgr- In this case 7^1/3 and so a return 
to a power law density pdf should occur for densities 
above (1 — 5) x lO'' cm~'^. However the grain, and 
hence gas, temperature is limited by the cosmic back- 
ground radiation to T w 3K, so it may be difficult to 
detect the decrease in gas temperature in this density 
range. 

If embedded protostellar sources are present, the 
grain temperature will probably exceed the gas tem- 
perature and gas-grain coUisional heating will dom- 
inate cosmic ray heating at large densities, as origi- 
nally emphasized by Goldreich and Kwan (1974). In 
this case it is easy to show that 7 should be signifi- 
cantly larger than unity at high densities. For exam- 
ple, at n = (10**, 10'^) cm~^, 7 = (1.4, 1.3), assuming 
c = 1.2 and Tgr S> T, and including the uncertain 
density dependence of the cooling function in eq. (|l^). 
These results are uncertain because a reliable treat- 



ment of the grain temperature requires a radiative 
transfer calculation. 

4. Extremely large densities. As pointed out by 
Lis & Goldsmith (1990) and others, the gas-grain cou- 
pling at very large densities is so strong that the gas 
temperature will simply follow the grain temperature, 
while the latter is controlled by the radiation field. So 
7 = 1 may obtain at very large densities when embed- 
ded protostars are present, if the ambient radiation 
field is not coupled to the local gas density. If positive 
feedback between density and star formation rate ex- 
ists, then the local ultraviolet energy density will scale 
similarly. The resulting grain temperature is only a 
weak function of energy density, scaling roughly with 
the 1/5 power. For example, if the local star forma- 
tion rate per unit volume scaled with the square of 
the density, d\ogT/d\ogn would be 0.4 and 7 would 
be 1.4. Accurate observational determinations of the 
gas temperature and density in high-density clouds 
could resolve the issue. 

A final remark is in order. It is well known that 
the Coriolis force (e.g., Binney & Tremaine 1987) and 
the magnetic field (e.g. Shu, Adams & Lizano 1987; 
Mouschovias 1987) may act against gravitational col- 
lapse. In fact, Vazquez- Semadeni et al. (1996) noted 
that inclusion of these effects restores a "higher-7" 
behavior to the numerical simulations, when consid- 
ering the ability of the flow to collapse gravitation- 
ally. However, the fact that in the present paper the 
density pdfs of simulations including both these ef- 
fects still exhibit power-law tails on their high-density 
sides suggests that the local production of large tur- 
bulent density fluctuations generated by the stellar 
energy input is not strongly inhibited by these pro- 
cesses. The effective compressibility of the flow prob- 
ably still has an equivalent polytropic exponent sig- 
nificantly smaller than unity. 

4. Applications 

4.1. Can the density pdf be used to derive 
the stellar IMF? 

Padoan, Nordlund, and Jones (1997, PNJ) have 
attempted to derive the mass spectrum of collapsing 
gas clouds from the density pdf of simulations. The 
procedure is as follows. If f{p) is the density pdf, then 
g{p) ^ pf{p) is the fraction of mass at a given density. 
Assuming that objects collapse if their mass exceeds 
the Jeans mass Mj, proportional to p^^^^, they pro- 
pose that the fraction of collapsing structures of mass 
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< M is proportional to J^°^ g{p)dp. The IMF is then 
the distribution of Jeans masses, obtained by differ- 
entiating this integral, giving an IMF proportional to 
g{p)dp/dMj. Since dp/dMj cx MJ^ , this gives an 
IMF which is an power law multiplied by the 

density pdf expressed in terms of Jeans mass, which 
they find from numerical simulations to be a lognor- 
mal. This latter factor is Gaussian in InM^, with a 
mean that depends on the temperature and turbulent 
velocity dispersion. This dependence causes the log- 
arithmic slope of the IMF to depend on these param- 
eters. In principle this procedure could be applied to 
any density pdf, such as those discussed above. Our 
concern here is the validity of using a one-point den- 
sity pdf to derive a distribution of masses, since mass 
as an attribute implies a coherent, contiguous "ob- 
ject" or structure, while the density pdf contains no 
spatial information of this nature. The problem is 
basically that not just any region of a given density 
can form a collapsing object, it must have a size large 
enough to contain a mass corresponding to at least 
the Jeans mass at that density. 

The problem can be illustrated using a specific ex- 
ample. Since Mj oc p~^l^ ^ the largest mass objects 
should form where the density is smallest. The low- 
est density regions are the "valleys" or "voids" in the 
density field. Imagine the density field as an array of 
pixels. If we consider a density value far into the low- 
density tail of the density pdf, a lognormal form of 
the pdf indicates that only a small fraction of pixels 
can have such small densities, and there is no a 'pri- 
ori reason to think that these pixels will be sufficiently 
contiguous that there really exists a coherent region of 
that density large enough to contain the correspond- 
ing Jeans mass (or any specified mass). In other 
words, in order to collapse, these low-density pixels 
must be spatially connected such that they form a 
very large low-density region. In general we expect 
the low density valleys to be scattered through the 
density field in some way. A Jeans' mass worth of 
matter centered at the position of a low-density pixel 
would be expected to contain pixels of a range of den- 
sities, some very large. This inconsistency occurs at 
any density-the low density case is just a severe ex- 
ample. Even at pixels where the density is largest, 
so the Jeans mass will be smallest, there is no reason 
to assume that nearby pixels out to mass Mj will be 
at that density]^ These same remarks apply to other 

^ A similar problem has been pointed out about Vazquez- 



instability criteria (since it is not at all clear that the 
thermal Jeans mass plays any role; e.g. Wiseman & 
Adams 1994; Simon 1997; Chappell & Scalo 1997). 

PNJ apparently recognized this problem, and state 
that 99% of the density peaks in their simulations 
contain mass exceeding the Jeans mass. However, 
our concern here is that, according to the formula- 
tion used by PNJ, the high-mass portion of the IMF, 
which is the region of most interest, is determined by 
low-density side of the density pdf. The problem then 
is that the low-density part of the pdf corresponds 
to voids in the density distribution, not peaks, and 
these low-density regions would have to be coherent 
over extremely large scales to contain a Jeans mass 
(which is largest at smallest densities). Furthermore, 
it is well-established observationally that stars form in 
the peaks of the density distribution (i.e., the clouds), 
not the voids. The implication that the highest-mass 
stars form in the lowest-density regions contradicts 
a large body of evidence on local star-forming com- 
plexes showing that high-mass stars form preferen- 
tially in the densest regions (e.g. compare Orion and 
Taurus, or the spatial distribution stars by luminosity 
within individual clusters in Orion and other regions 
that contain massive stars). 

Yet another concern is the statement itself by PNJ 
that the vast majority of their density peaks is Jeans- 
unstable, since it has been suggested on both obser- 
vational (e.g., Falgarone, Puget & Perault 1992; Mag- 
nani, LaRosa & Shore 1993) and numerical (Vazquez- 
Semadeni et al. 1997) grounds that a large fraction of 
the turbulence-induced transient density fluctuations 
are not gravitationally bound. Besides, PNJ's simu- 
lations are not reported to be self-gravitating, so the 
statement that the density peaks are unstable appears 
inconsistent. 

We conclude that the one-point density pdf can- 
not be used to derive a mass function of "clouds" 
or protostars: spatial statistics are necessary. One 
might consider the two-point density pdf, which gives 
the probability that pairs of points separated by r 
have densities pi and p2- A second-order moment 
of this two-point probability distribution, the corre- 

Semadeni's (1994) criterion for the development of hierarchical 
structure based on the density pdf (A. NouUez, private commu- 
nication). Nevertheless, in that work, an additional constraint 
was required from the density pdf, namely, that it exhibited 
self-similarity as smaller and smaller regions are considered. 
This requirement may possibly provide the necessary spatial 
information. 
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lation function, reduces this information greatly (in 
the same way that the variance reduces the one-point 
pdf to a single number) . But even the 2-point pdf is 
insufficient to estimate the mass spectrum because we 
need to know the probability that there exists a re- 
gion with enough contiguous pixels in a given density 
range to comprise a mass sufficient for instability. At 
the least, a joint probability distribution of densities 
and coherence sizes would be needed, but construc- 
tion of this joint distribution is equivalent to directly 
constructing the mass spectrum from the simulations 
or observational data. The plausibility of our conclu- 
sion can also be seen by considering the process of 
estimating a mass spectrum for interstellar "clouds" 
from observational data. Within an arbitrary density 
field (or a three dimensional field with radial velocity 
as one axis, as for spectral line data) the identification 
of entities defined in some operational manner and the 
measurement of their attributes is a complicated and 
subjective procedure (e.g. Houlahan and Scalo 1992, 
Williams et al. 1994), while the estimation of the col- 
umn density pdf by histogram construction for the 
same field is simple by comparison. This difference 
in computational effort is a reflection of the essential 
role of spatial information in the former procedure 
and the absence of such information in the latter. 

It is interesting to note that the procedure pro- 
posed to derive the mass spectrum of bound conden- 
sations by PNJ is very similar in approach to the 
derivation of the mass spectrum of bound objects in 
a gaussian density field (in the linear regime) used 
in considerations of galaxy formation (e.g. Press & 
Schecter 1974, Peacock & Heavens 1990). There the 
threshold density is taken as the critical density con- 
trast needed at some initial time so that the contrast 
will be about unity at a later time. But in the cosmo- 
logical case, the derivation assumes that the density 
field is succesively smoothed by window functions of 
scale R, which would correspond to a mass M pro- 
portional to the mean density in the window. This 
smoothing takes care of the problem discussed here, 
of having enough extent at each critical density to 
contain mass M, but the answer, even for a gaussian 
random field, depends on the power spectrum and the 
rather arbitrary choice of window function. In the 
case of simulations the problem is more difficult. In 
particular, we suspect that because the density field is 
non-gaussian, the answer will depend on higher-point 
correlations, as we speculated earlier. 

In any case, using simulations alone, we think that 



the only way to determine the mass spectrum of co- 
herent condensations, whether collapsing under the 
Jeans criterion or not, is to directly identify "clouds" 
in the simulation and calculate their frequency distri- 
bution (e.g. Vazquez-Semadeni et al. 1997; Chappell 
and Scalo 1997). There is no "shortcut" to the IMF 
through the one-point density pdf. 

Besides the conceptual problem of deriving an IMF 
from the one-point density pdf without spatial in- 
formation, the IMF derived by PNJ docs not give 
a good representation of the observed IMF in either 
field stars, open clusters, OB associations, or nearby 
galaxies. The PNJ IMF decreases with increasing 
mass with a power law index that decreases (becomes 
more negative) with increasing mass. PNJ recognized 
that the function does not match the observed IMF 
at any single gas temperature and velocity dispersion, 
and so attempted to remedy the situation by assum- 
ing a T^^ distribution of gas temperatures between 5 
and 40 K and showed that such a temperature distri- 
bution could approximately match the Miller & Scalo 
(1979) field star IMF. The problem is that the Miller- 
Scalo IMF at higher masses has long been superseded 
by a large body of observations (see Scalo 1986 for 
a summary of work before 1985; also Rana 1991 for 
field stars, Massey et al. 1995a, b for massive stars in 
the Galaxy and Magellanic Clouds; see Scalo 1997 for 
a recent review) that suggest that for masses above 
about 1.5 Mq, the IMF is not lognormal, but is best 
fit by a power law of slope around —1.7 ± 0.2 with 
(probably) a fiatter high-mass tail. Such a form would 
be difficult to produce from the PNJ formalism with- 
out a very contrived mixture of gas temperatures and 
velocity distributions.^ 

For the above reasons, we must conclude that 
the application of the PNJ IMF to globular clus- 
ters (Padoan, Jimenez, and Jones 1997), low surface 

*For comparison, if the density pdf were a power law of the form 
f{p) ~ , with a sharp cutoff at small densities, the IMF 
for a single-temperature Jeans mass would be proportional to 
M^^~^, or M~^'^ for 9 = 1.7 as found here, in good agreement 
with observations. The flattening at large masses could be pro- 
vided by the flattening of /(p) at densities near the peak. The 
turnover of the IMF at small masses could be ascribed to "tur- 
bulent viscosity" steepening the density pdf at large densities 
(which requires 9 > 2.5 in this regime if it were a power law). 
By assuming a parameterized distribution of temperatures, we 
could match essentially any desired form of the IMF. However 
we emphasize that this procedure is unwarranted, and only give 
the example as an illustration of the ease with which such a pa- 
rameterized IMF model could produce a match to any set of 
observations. 
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brightness galaxies (Padoan, Jimenez, and Antonuccio- 
Delogu 1997), and other galaxy evolution problems 
(Chiosi et al. 1997) are without basis. The galaxy 
evolution models of Chiosi et al. (1997) do represent 
a major advance in such modelling, in that the cou- 
pling between the IMF and the energy equation is 
explicitly treated. However we interpret their results 
as showing that an IMF that favors higher mass stars 
for some combination of lower mean gas density and 
higher temperature and/or velocity dispersion can ac- 
count for a number of observed chemical and pho- 
tometric constraints for elliptical galaxies. We only 
point out that there is as yet no physically consistent 
model for such an IMF and that other models besides 
the PNJ model could be constructed to give similar 
dependences. 

We emphasize that the fact that a theoretical for- 
mulation of a problem has a large enough number 
of adjustable parameters to interpret a number of 
diverse phenomena is not a measure of its physical 
validity. A theoretical IMF that can be applied to 
these and other problems must be able to give a nat- 
ural account of contemporary studies of the local field 
star and cluster IMF, and, besides, cannot be derived 
solely from the density pdf. A brief summary of re- 
cent theoretical models for the IMF is given in Scalo 
(1997). 

4.2. Applications to Column Density Statis- 
tics 

Our critique of the PNJ IMF does not extend to 
their interpretation of the increase of extinction stan- 
dard deviation with average extinction found by Lada 
et al. (1994), given by Padoan, Jones, & Nordlund 
(1997). We agree that this increase constrains the 
density pdf and probably requires an intermittent tail 
(not necessarily lognormal, though) and a power spec- 
trum of the density field P{k) oc fc^, with x around 
—2 to —3. However, whether or not simulations with 
different physical assumptions can give the requisite 
intermittency and power spectrum remains an open 
question. We have examined the power spectrum of a 
2-dimensional MHD model (fig. 12) and find a small 
intermediate range of wavenumbers over which the 
power spectrum is a power law with index around 
—2.4, close to the value of —2.6 reported by Padoan, 
Jones, and Nordlund (1997) in their simulations, al- 
though a much flatter and just as extensive power law 
occurs at smaller wavenumbers in the present simu- 
lations. A more direct test of the simulations would 



be an estimate of the pdf of column densities. Such a 
comparison, using IRAS 100 fim column densities for 
low-mass star regions and ^^CO data for high-mass 
star regions, is postponed to a separate paper (Scalo, 
Chappell, Miesch, and Vazquez-Semadeni, in prepa- 
ration) . 

Note that in our simulations, the wavenumber 
range with a —2.6 "slope" may actually constitute 
the range at which the mass diffusion is active, pro- 
ducing an extended exponential decay which may be 
confused with a power-law. This effect could also be 
at work in the simulations of Padoan, Jones, & Nord- 
lund(1997), due to the numerical diffusion. Note, 
however, that this cannot be at the origin of the 
power-law tails we observe in the density pdfs from 
our simulations, since the effect of the diffusion term 
is to smooth oiit fluctuations, while the power-law tails 
in the pdfs are enhancements over the lognormal al- 
ternative. 

4.3. Size distribution of wind-driven shells 

Shells driven by H II regions, protostellar winds, 
supernova remnants (SNRs), and superbubbles are 
believed to play an important role in shaping inter- 
stellar structure over a range of size scales, and the ve- 
locity and size distributions of such shells are of inter- 
est in a variety of astrophysical contexts, such as the 
porosity of the ISM (Oey & Clarke 1997) and models 
for the IMF (Silk 1995). The expansion law for shells 
depends on the ambient interstellar density, so a dis- 
tribution of ambient interstellar densities could affect 
the predicted size and velocity distributions. The dis- 
tribution of sizes of SNRs and superbubbles for a dis- 
tribution of source luminosities has been treated in 
detail by Oey & Clarke (1997), who assumed a con- 
stant ambient density. Here we illustrate the role of 
a density pdf by assuming a constant luminosity. For 
simplicity, we neglect the effect of shell stalling which 
occurs when the driving pressure matches the ambi- 
ent pressure. We are mostly interested in effects that 
occur in molecular clouds subjected to internal proto- 
stellar winds. In this case the appropriate expansion 
law is probably that for a non-adiabatic momentum- 
conserving shell (see the conditions derived in Nor- 
man & Silk 1980), which can be expressed as 



(15) 



where L is the source luminosity, which we assume is 
constant for illustrative purposes, and R the radius of 
the shells. 
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For a constant star formation rate £>, the size 
distribution for a given ambient density is given by 
N{R, p)dR = Bdt, so, integrating over the pdf of am- 
bient densities, 



NiR) 



fip){dR/dt)-^dp. 



(16) 



If the ambient densities were the same for all 

sources (/(p) a delta function), then, since eq. (15) 
gives ipR/di)^^ ~ i?, N{R) would scale as R, reflect- 
ing the larger number of shells at small velocities. For 
the case of a distribution of ambient densities, the in- 
tegration limits require some care. We assume that 
the lower and upper limits of ambient densities are 
Pmin and Prnax, and that all the soTirces are only ac- 
tive for a time t^- Then for radii R less than some 
small R — Ri, all ambient densities, even up to the 
maximum density, allow expansion to radius i? in a 
time tg. But for larger radii, shells expanding into 
high density material will not reach that size in time 
tg, so the limits of integration are pmin to p{R, te), 
where p{R, tg) is the ambient density for which a shell 
will reach size R within time tg- So 

fPiR, *e) 

N{R) ~ / f{p)p'/^Rdp. (17) 

J Pmin 

Taking f{p) ^ gives 



N{R)^ p{R,tg) 



-e+3/2 



-0+3/21 „ 
Pmin 



(18) 



For > 3/2 (steeply declining pdf), most of the ambi- 
ent material is at small p, and the second term domi- 
nates. The size distribution is then unaffected by the 
density pdf, N{R) ~ R, because the pdf is in effect a 
delta hmction. 

But if 9 < 3/2 (flatter pdf), the first term domi- 
nates. Since p{R,tg) for a fixed tg, the result- 
ing size distribution is 



N{R) ~ R 



(19) 



For 6 < 3/2, 46—5 < 1, so N{R) increases less rapidly 
than R; shells at small R become significant because 
more of the ambient material is at large densities. 
The same derivation for an adiabatic pressure-driven 
shell (Weaver et al. 1977) gives N{R) ~ ijse-s/s ^j. 
e > 2/3. 

The situation is actually much more complicated 
for a number of reasons. The density pdfs found in 



the simulations are strongly decreasing functions of 
p for large densities, but increasing, and extremely 
uncertain, functions of p at small densities. Numeri- 
cal integration of eq. (17) for a double power- law or 
lognormal density pdf would be a useful exercise, but 
our purpose here is only to illustrate the potential 
effect. Secondly, a distribution of source luminosi- 
ties should be included (i.e. a generalization of Oey 
fc Clarke's 1997 result to include the density pdf). 
Third, as noted above, shell interactions are probably 
important (see also Norman & Silk 1980), and for this 
process the effect of a distribution of ambient densi- 
ties woiild be to affect the column density and veloc- 
ity distributions of the interacting shells. Finally, we 
note that the present derivation suffers from the same 
lack of spatial information that we claimed precludes 
a derivation of the mass spectrum. If the ambient 
density fluctuations were all on scales smaller than 
any shell size of interest, then the effect of the den- 
sity distribution would only be to corrugate the shell 
in some irregular manner, as different parts of a given 
shell expand into different ambient densities. On the 
other hand, some coherent density fluctuations in the 
ISM have scales larger than the shell sizes of interest. 
A full examination of these problems is beyond the 
scope of the present paper. 

4.4. Comparison with cosmological simula- 
tions 

It is of interest to compare our results with cosmo- 
logical simulations of the evolution of density fluctua- 
tion in the universe. Such calculations usually contain 
a minority of the matter in baryonic form, with the 
rest either in "cold" or "hot" dissipationless particles 
meant to represent varying dark matter candidates, 
or some combination of the two. The dissipationless 
component is equivalent to a self-gravitating Burgers 
flow (no pressure, 7 = 0) with viscosity arising only 
from numerical diffusion, while the coupled baryonic 
component follows the Navier-Stokes equation with 
self-gravity. 

Most of the work in this held has been aimed at 
using the density pdf as a discriminator of parame- 
ters related to initial conditions, such as the index of 
the initial power spectrum or deviations from initial 
Gaussian statistics. Several papers, apparently begin- 
ning with Hamilton (1985), have proposed that the 
density pdf in the not-too-strongly nonlinear regime 
is lognormal (see Weinberg and Cole 1992, Kofman 
et al. 1994, Bernardeau and Kofman 1995; for other 
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references see Protogeros and Scherrer 1997). How- 
ever Bernardeau and Kofman (1995) argued that it is 
likely that the lognormal distribution only obtains for 
an initial power spectrum with index not too different 
from —1 (as in cold dark matter scenarios) and for a 
variance of the (smoothed) density field not too large 
(less than around unity). For conditions outside of 
this range the density pdf appears to have a power- 
law, not lognormal, high-density tail (Bouchet and 
Hcrnquist 1990, Melott et al. 1997; these power laws 
are to be distinguished from the steeper power laws 
that result from caustics in the Zeldovich approxima- 
tion with no smoothing). Although in general such 
simulations are not isothermal and in some cases refer 
to pressureless dark matter, the apparent transition 
from lognormal at small density variance to power law 
at larger density variance may be understood as the 
analogue of the Mach number dependence described 
in the present paper, the range of density values being 
just too small for the power-law range to develop. 

A recent discussion of a cold-plus-hot dark mat- 
ter simulation that presents the density pdf has been 
given by Klypin, Nolthenius, and Primack (1997). 
These simulations use a particle-mesh code with a 
force mesh resolution of 512'^ in three dimensions. 
The resulting function pf{p) (their fig. 16) shows a 
clear power-law over more than two orders of magni- 
tude in density, with a power-law index of about —1.2, 
and a dropoff at the highest densities due to finite 
resolution. There is certainly no hint of a lognormal 
form. The result supports our contention of a power 
law density pdf driven by advection, and shows that 
power-law behavior occurs even in three dimensions. 
The power law index is steeper than found here by 
about 0.2 to 0.6, but it is uncertain whether this is a 
result of the different dimensionality, the fact that the 
self-gravity is mainly supplied by the dissipationless 
component, or some other effect. 

5. Conclusions 

In this paper we have discussed various theoretical 
aspects of the density probability distribution func- 
tion, or pdf, of the ISM. We first presented evidence 
that the pdf appears to be a power-law at high densi- 
ties in a number of two-dimensional numerical simu- 
lations of the ISM at intermediate scales, in which the 
various intervening physical agents are progressively 
removed until only the nonlinear advection operator 
of the momentum equation is left. Additionally, in all 



cases the morphology is extremely filamentary. This 
result has two main implications. First, it suggests 
that both the functional form of the pdf and the mor- 
phology of the density field are the signature of the 
advection operator. 

Second, this result is in contrast with previous find- 
ings that the pdf is lognormal (Vazquez-Semadeni 
1994; Padoan et al 1997 a, b, c). To resolve the dis- 
crepancy, we turned to the one-dimensional simula- 
tions of PVS, which suggest that the lognormal pdf 
is realized only in the isothermal case (7 = 1), while 
at smaller values of 7 a power law appears. More- 
over, PVS report that the deviation from a lognor- 
mal is larger at larger Mach numbers. These results 
are consistent with the fact that the simulations of 
Vazquez-Semadeni (1994) and Padoan et al. (1997a, 
b, c) were isothermal, while those presented here in 
all cases have small or zero values of 7. 

The dependence on the polytropic index 7 prompted 
an investigation on what are its expected values in the 
actual ISM, as determined by the equilibrium between 
heating and cooling rates, which in turn depends on 
the density and temperature of the medium. From 
inspection of published heating and cooling rates, we 
expect that at densities below ~ 10^ cm~^ values of 
7 significantly smaller than unity should occur, while 
at larger densities, values near unity should appear. 
The switch-over occurs because of the effect of radia- 
tive trapping in decreasing the density dependence of 
the molecular cooling rate at densities above ~ 10'^ 
cm~^. At densities larger than around 10^ cm~^ gas- 
grain coUisional cooling or heating should dominate, 
forcing 7 back away from unity. Thus, lognormal pdfs 
are only expected in the range 10^ to (1 — 5) x 10"^ 
cm~^, with power law pdfs outside of this range of 
densities. However, at extremely large densities the 
gas-grain coupling may be so efficient that the gas 
temperature will simply follow the grain temperature 
which is controlled by the ambient UV radiation field. 
To the extent that the radiation field energy density is 
independent of the gas density, the temperature will 
then be constant, so a return to a lognormal density 
pdf may occur at such large densities. A star forma- 
tion rate-density coupling would increase 7 somewhat 
above unity through the effect of the young massive 
stars on the grain temperature. 

A nmnber of possible applications were discussed. 
Most importantly, we discussed the feasibility of de- 
riving the mass spectrum of bound condensations or 
the stellar IMF from knowledge of the density pdf 
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alone, as recently done by PNJ. We conclude that 
this is not possible, because the pdf does not contain 
spatial information, while the gravitational instabil- 
ity criterion requires simultaneous knowledge of mass 
and size information. Second, we briefly discussed the 
applicability to extinction statistics (Padoan, Jones & 
Nordlund 1997). Third, we analyzed the effect of an 
inhomogenous density field on the size distribution of 
wind-driven shells in the ISM. We find that in a sim- 
ple example case with f{p) ^ , the resulting dis- 
tribution increases more slowly with shell size than in 
the case of a uniform density background if 6 < 3/2. 
However, we warn that the real situation is likely to 
be much more complicated. 

Finally, we compared our results with those from 
a number of cosmological simulations, in which both 
lognormal and power-law forms of the pdf have been 
reported. We speculate that the appearance of log- 
normal pdfs for non-isothernial situations may be a 
consequence of small density variances, analogous to 
the Mach number effect described here. 

A central result of our work is the conclusion that 
it is primarily the nonlinear advection operator oper- 
ating on a source of compressible motions (whether it 
be self-gravity or nonlinear waves or cooling or some- 
thing else) that is responsible for both the general 
filamentary morphology (which would include sheets 
in 3D) and for the power law portion of the density 
pdf. This result is consistent with the presumed ten- 
dency of the advection term to induce multiplicative 
processes which lead to power law scaling, as it does 
in the incompressible case. If this is correct, then the 
simulations studied here should lead to investigations 
of analytical models for the behavior of the advection 
operator in a highly compressible medium. For that 
purpose our simulations provide a constraint in the 
form of the density pdf that must be explained both 
when 7 is and is not close to unity, as well as a con- 
straint on the morphology. As an example, we point 
out that Elmegreen's (1997) construction of a fractal 
by nesting can be shown to produce a density pdf that 
has the form f(p) ^ P^^, independent of fractal di- 
mension, which is not much different from the power 
laws found here. Although Elmegreen's model is a 
static construction with no reference to underlying 
physics such as advection, this does suggest that an 
advection model that is related to hierarchical nesting 
can explain the power-law density pdf. 

Alternatively, the physics might be multiplicative 
but not involve nesting. For example, the action of 



advection could simply be to steepen fluctuations into 
"shells" , or "filaments" or "spikes" (depending on di- 
mensionality), as occurs in the Burgers equation. Af- 
ter one crossing time the evolution is dominated by 
the merging of these entities (which we might refer 
to as "blobs"), a process that is described by the so- 
called "coalescence equation" . From much previous 
work we know that the solution of this kinetic equa- 
tion, when forced at small scales, is a power law mass 
spectrum with index around —1.5 to —2. If the blobs 
were spikes or filaments or shells whose thicknesses 
were relatively constant or at least independent of 
their mass, the density pdf would be a power law in 
the same range. Obvious variations and elaborations 
on this theme suggest themselves. 

Our point is that it is relatively easy to think of 
simple models that might explain a power law den- 
sity pdf. However much more work is needed to un- 
derstand the actual physics of advection acting on a 
source of compressibility, and how that physics gives 
rise to the observed and simulated morphology and 
statistical properties. Given the longstanding diflacul- 
ties with the analogous questions in the incompress- 
ible case (in which vorticity dynamics plays the key 
role), we do not expect oversimplified models to pro- 
vide the necessary answers. 
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Fig. 1. — Logarithmic grayscale image of the density 
field of run MHD800 at t = 1.2io, where to = 8.2 x lO'^ 
yr is the sound crossing time. The simulation repre- 
sents a region of size 1 kpc, has an initial temperature 
of 10^ K, and has an initially turbulent velocity field 
with an rms velocity fluctuation of 11.7 km s~^. The 
maximum and minimum values of the density at the 
time shown are Pmax = 526.3 and Pmin = 2.42 x 10~^. 



This 2-column preprint was prepared with the AAS lAT^^X 
macros v4.0. 
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Fig. 2. — Histograms of log p for run MHD800 at a) 
t = 0.9io (left) and b) t = 1.2io (right). The localized 
stellar heating was turned off at i = QMo in a) and at 
t = 1.16fo in b). The density field at the latter time is 
shown in fig. 1. In a) a clear power law is seen in the 
range —0.6 < logp < 0.9, with slope ^ —0.73, while 
in b) a slightly less well-defined power law is seen for 
< logp < 1.5, with slope ^ —0.63. 



Fig. 4. — Histograms of logp for run HD512 at times 
a) t = 0.4^0 (left) and b) t = I Alto (right), respec- 
tively shown in figs. 3a and b. In a) two power-laws 
with slopes —0.65 and —1.9 can be seen, although 
the plot can also be interpreted as a single power law 
with a bump at logp 0.3. In b) a single power law 
with slope ~ —1.0 is observed, with two bumps, one 
at logp ~ 0.7 and the other at logp ~ 2.7. These 
seem to be respectively due to the onset of star for- 
mation at p = 30 (see text) in the first case, and to 
the collapsed region, whose peak density is p = 703. 






Fig. 3. — Density fields for run HD512 at times a) 

t = 0.4io (left) and b) t = 1.47to (right). This is a 
non-magnetic run otherwise similar to run MHD800 
(except for the resolution). In a) a network of fila- 
ments is observed, while in b) material in the upper 
half of the run has collapsed gravitationally. Star for- 
mation remained on at all times in this run, although 
it was not able to avoid the gravitational collapse ob- 
served in b). 



Fig. 5. Spatial distribution of the square root of 
the density field at times 0.1 (left) and 10 (right) 
initial crossing times for a 256^ quasi-Burgers decay 
run. The simulation was initialized with a constant 
density and a velocity field with energy spectrum 
E{k) a fc*exp(— A;^/fco) with ko = 8/L, where L = 
size of simulation box. 
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Fig. 6. Time development (top to bottom) of pf{p) 
for two 256^ quasi-Burgers decay runs. Initial energy 
spectra are proportional to (left) and fc" (right). 
Times, from top to bottom, are 0.006, 0.07, 1.1, and 
17 initial crossing times. The solid lines are reference 
lognormal pdfs with peaks at p = 1. 
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Fig. 7. Late-time distributions pf{p) for three 256^ 
quasi-Burgers decay simulations initially excited at a 
single wavenumber fco = 64 (top), 16 (middle), and 
4 (bottom). Solid lines are reference lognormal pdfs 
with peaks at p = 1. 
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Fig. 8. — Distribution pf{p) for a 512^ quasi-Burgers 
decay simulation at times 0.1 (top) and 2 (bottom) 
initial crossing times. Solid lines are reference lognor- 
mal pdfs with peaks at p = 1. 
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Fig. 9. — Density fields for two onc-dimcnsional runs 
with M = 3 (see eq. [12]), with a resolution of 2048 
grid points, a) 7 = 1 (top) and b) 7 = 0.3 (bottom). 
These runs all start from rest, but are forced with 
white noise at the large scales (see text) and with a 
correlation time tcorr = 4.77 x 10~*to- The run with 
7 = 0.3 shows taller, narrower density peaks than 
that with 7 = 1. 
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Fig. 10. — Histograms of log p for the ID simulations 
shown in fig. ^ with M = 3. The runs are sam- 
pled at intervals At — 1.59 x lO^^to over a total time 
t = 23.9fo- The histograms thus include 1500 x 2048 
data points. The dotted lines show lognormal fits to 
the curves, a) 7 = 1 (top). The fit by a lognormal 
is excellent, b) 7 = 0.3 (bottom). The histogram 
differs drastically from a lognormal, and instead a 
clear power-law is seen at high densities, in the range 
0.3 < logp < 1.3. 
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Fig. 11. — Histogram of logp for a ID simulation 
similar to that shown in fig. |^a, except with M = 
1.68. The deviation from a lognormal is seen to be 
less pronounced than in the case M = 3. Also, the 
power-law at high densities is less well-developed. 
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Fig. 12. — Density power spectrum for run MHD800 
at time t = 0.9to- Two power laws are seen, one in the 
range 0.6 < logfcl.5 (fc is the wavenumber) with slope 
~ —0.9, and the other in the range 1.5 < log A: < 2.3, 
with slope —2.4. The region with the steeper slope 
may actually not be a true power law, but rather a 
slow exponential decay due to the mass diffusion term 
in the continuity equation. 
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